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Induced Compton scattering (ICS) is an interaction between intense electro-magnetic 
radiations and plasmas, where ICS transfers the energy from photons to plasmas. 
Although ICS is important for laser plasma interactions in laboratory experiments and 
for radio emission from pulsars propagating in pulsar wind plasmas, the detail of photon 
cooling process has not been understood. The problem is that, when ICS dominates, 
evolution of photon spectra is described as a nonlinear convection equation, which makes 
photon spectra to be multi-valued. Here, we propose a new approach to treat evolution of 
photon spectra affected by ICS. Starting from the higher-order Kompaneets equation, we 
find a new equation that resolves the unphysical behavior of photon spectra. In addition, 
we find the steady-state analytic solution, which is linearly stable. We also successfully 
simulate the evolution of photon spectra without artihcial viscosity. We find that pho¬ 
tons rapidly lose their energy by ICS with continuously forming solitary structures in 
frequency-space. The solitary structures have the logarithmically same width character¬ 
ized by an electron temperature. The energy transfer from photons to plasma is more 
effective for broader spectrum of photons such as expected in astrophysical situations. 
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1. INTRODUCTION 

Nonlinear interactions between strong electromagnetic waves and plasmas have been studied 
in the context of laser plasma physics and astrophysical phenomena. Depending on intensi¬ 
ties of radiations and conditions of plasmas, plenty of nonlinear interactions exist and are 
studied in different approaches. For example, quantum electrodynamics predicts vacuum- 
polarization effects, where photons interact with virtual electron-positron pairs (c.f. Ref [1]). 
Parametric processes in plasmas, such as induced scattering processes, filamentation insta¬ 
bility and so on, have been studied in classical electrodynamics with the two-fluid description 
of plasmas (c.f. Ref. [2]). There are also studies of parametric processes in the semi-classical 
formulation, which is used in this paper (c.f. Refs. [3-5]). Here, we focus on induced Comp¬ 
ton scattering (ICS), which is a parametric instability between photons and each electron, 
rather than plasmon. The condition for ICS to be dominant process is when both the central 
frequency and the spectral width of strong electromagnetic radiations are greater than the 
Langmuir plasma frequency [6]. ICS cools photons toward the Bose-Einstein condensation, 
while plasmas are heated by this process (c.f. Ref. [7]). 
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ICS has been studied in the context of laser plasma interaction both theoretically [8-10] and 
experimentally [11-14]. Most of them studied the process with the two-fluid approximation, 
which accounts only for the linear regime of the instability. On the other hand, ICS has also 
been studied for the scattering process around high intensity astrophysical sources with the 
semi-classical formulation, which can partly treat the nonlinear regime (e.g. Refs. [15-18]). 
However, they constrained plasma conditions around some astrophysical objects considering 
only the scattering optical depth to ICS. It is known that there is a difficulty to treat the 
nonlinear regime of ICS even with the semi-classical formulation as described in this paper. 

In the semi-classical formulation, a relaxation process of isotropic photons interacting 
with a rarefied Maxwellian plasma by Compton scattering is described by the Kompaneets 
equation [31]. The Kompaneets equation is the Fokker-Plank equation describing evolution 
of a photon spectrum including the ICS term, which comes from the Boson nature of photons 
and, which is quadratic in the intensity of the radiation. Therefore, for high intensity radi¬ 
ations, this quadratic term plays a dominant role and the Kompaneets equation is reduced 
to the nonlinear convection equation (Eq. (22)). It is known that the nonlinear convection 
equation has the implicit solution, which will be multi-valued after a finite time starting 
from certain initial conditions. Such a spectrum is unphysical and should be bent to be 
single-valued, i.e., there should be an appropriate physical process which we have ignored 
(e.g. Ref. [19]). 

Different approaches have been adopted to resolve the unphysical behavior of a photon 
spectrum when ICS dominates. Peyraud (1968) [20-22] considered heuristically to add a dis¬ 
persive term (the third derivative term) to the nonlinear convection equation and made the 
equation Korteweg-de Vries-type, i.e., soliton formation in frequency-space. Reinish (1976) 
[23, 24] considered to recover the diffusive term (the second derivative term), which is orig¬ 
inally included in the Kompaneets equation, and made the equation Burgers-type, which 
allows formation of a shock wave structure like hydrodynamics. However, because the orig¬ 
inal diffusive term is linear in the intensity of radiations, this will not be applicable when 
the intensity of radiations becomes higher and higher. On the other hand, Zel’dovich and 
colleagues [25, 26] adopted the integral form of the equation, i.e., the Boltzmann equation 
(Eq. (1)). They predicted formation of solitary structures in a radiation spectrum. It is worth 
distinguishing which spectral behaviors appear in a given condition between photons and 
electrons. There are also numerical studies of this problem by Montes (1977, 1979) [27, 28] 
and Coppi et al. (1993) [29]. We compare their results with ours in Section 4.3. 

In this paper, we propose a new approach to study evolution of a radiation spectrum when 
ICS dominates. In Section 2, we extend the Kompaneets equation to the higher-order terms. 
This is a preparation to obtain a new equation and its derivation helps understanding what 
is the resolution of the problem in the nonlinear convection equation. In Section 3, we find 
a new equation that describes evolution of a radiation spectrum for intense radiations. We 
obtain the analytic solution in steady state and also discuss its linear stability. In Section 4, 
we show some examples of numerical studies of the new equation. We give interpretations 
of our results and comparison with past studies. Section 5 is devoted to summary. 

2. Higher-Order Kompaneets Equation 

Here, we reconsider a kinetic equation that describes the evolution of the photon occupation 
number n{k) by Compton scattering. Considering a rarefied plasma, we neglect emission 
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and absorption of photons by plasmas, i.e., the photon number is conserved in this system 
(c.f. Ref. [7]). Effects of the background magnetic field and polarizations of photons are 
also neglected. Below, a wavenumber k is in unit of the inverse of the electron Compton 
wavelength Ae = h/rrieC and a momentum p of an electron is in unit of m^c. We start from 
the Boltzmann equation for photons including Compton scattering by free electrons of a 
density ng. We use a normalized momentum distribution of electrons /(p), i.e., f d^pf{p) = 
1. Assuming the spatial homogeneity of the system, we have 


dn{k) 

~dr 


CUe / (fpfip) 


d?ki 


Din{ki){l + n{k)) 


Dn{k){l + n(fci)) 


da 


dkidQi 


/ /ci \ ^ da 
\ A; y dkdCl 


( 1 ) 


where the terms 1 + ra represent spontaneous and induced scattering terms, respectively. 
D = 1 — (3 ■ ft and Di = 1 — P ■ fii represent the factor due to relative velocities between 
interacting photons and electrons, and {ki/k)‘^ comes from the difference in phase space 
volumes between d^ki and d^k, where Cl = k/k, fti = ki/ki, and |/3| = (1 + 
respectively. The differential scattering cross section for Compton scattering (the Klein- 
Nishina cross section) from an initial state fc, to a final state fcf of a photon is (e.g. Ref. 
[30]) 
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where quantities with ~ represent those in the electron-rest frame, 7 = ( 1 -|-is the 
Lorentz factor of an electron, p = fli ■ Clf is the cosine of the angle between fci and fcf, 6 is 
the Dirac’s delta function, and cjrp is the Thomson cross section. We have used the following 
Lorentz transformation laws 


ki^f = 7A,fA:i,f, dftf 
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and {D\kj)/{Diki) + {Dikp/{D\k\) — 1 = 1 -|- kiki{\ — p^)/D^Df) with the use of 5[ki — 
'jDiki/i'jDi + ki{l - A)]- 

We execute the integral with respect to ki by using the 6 function in Eq. (2) and intro¬ 
duce notations y = ©ngO-rpCt, 0 = k^Tpi/nieC^, k± = ^Dk/{'yDi =p A:(l — p)), and n{k±) = 
n{k±,fti), where ks is the Boltzmann’s constant and Tp\ is an electron temperature. Then, 
Eq. (1) is expressed as 
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AS* and AI correspond to contributions from spontaneous and induced scattering, respec¬ 
tively. We find that net contributions to scattering are differences between an upper state 
and the incident state k for spontaneous scattering, and between and a lower state 
fc_ for induced scattering. 

Now, we assume that electrons and photons are non-relativistic, i.e., p 1 and A: <C 1. 
The lowest order in p and k gives Thomson scattering where k± ~ k, n{k±) ~ n{k, fii), and 
R^ ^ 3(1 -|- /i^)/167r. The result is 


ASo = {l + p^){n{k,ni)-n{k,il)), (8) 

AIq = {l + p^)n{k,fl){n{k,fli)-n{k,fti)) = 0, (9) 


where AS'(^) and A/(^) represent terms on an order ^ for AS" and AI, respectively. When the 
photon distribution is isotropic, we obtain A5o = 0. On the other hand, AIq = 0 is always 
satished because ICS requires the frequency shift in essence. 

On the hrst order with an assumption of isotropy, we obtain the well-known Kompaneets 
Equation. The right-hand side of Eq. (4) becomes 

f A'Kp^f{p)dp ^ ^ ^dflr 
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where we have used / Airp^f{p)dp = 30, i.e., the Maxwell-Boltzmann distribution. Though 
the result is simple, the derivation is fairly long. We just note the expansion of k± on the 
first order in k and 

k 

where a = Vt 


1 + p{ai — a) +p'^ai{ai — a) ± k{l — p), 


( 11 ) 


p • fl and ai = flp • fii. Odd orders in p will vanish because of isotropy 
(f dflidnpa^"^~^af~^ = 0 with natural numbers m and 1) and an order of p^ corresponds 
to the first order of the plasma energy, i.e., the temperature 0. Introducing normalized 
frequency x = k/Q, Eqs. (4) and (10) give the Kompaneets equation [31] 


dn{x) 

dy 




( 12 ) 


In addition, multiplying Equation (12) by x^ and integrating over dx, we obtain evolution 
of the photon energy density (cph = f x^n(x)dx) as (c.f. Ref. [32]) 

deph(y) 


dy 


= — j x'^n{x)dx — J x‘^n^{x)dx + 4:€pii. 


(13) 


The hrst term, which comes from ASy^y of Eq. (10), means energy loss of photons by spon¬ 
taneous Compton scattering. Corresponding induced term (A/^^)) is the second term, which 
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also implies energy loss of photons. On the other hand, the last term corresponding to A5(p2) 
expresses energy gain by inverse Compton scattering, even though this term comes from the 
elastic approximation (the Thomson limit, zeroth order in k). It should be noted that there 
is no induced term associated with inverse Compton scattering, i.e., = 0. As is well- 

known, this Kompaneets equation is not enough to follow evolution of the photon spectrum 
for 1. 

Now, we proceed to the next order. We have A5(fc2), AS'(fcp 2 ), and A5(p4) for spontaneous 
scattering but only a cross term exists for induced scattering, because the numerator 

of Eq. (6) vanishes for even orders in each k and p. The expansion of k± on the second order 
in k and (0) is 
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Note that only the terms oc k have double sign on the right-hand side of Eq. (14) and this 
is the reason for AI(fc 2 ) = A/(p4) = 0, i.e., odd orders in k is required for induced scattering. 
The integrals are carried out as 
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for induced scattering. We also need the relativistic correction to the Maxwell-Boltzmann 
distribution f{p) = e“'’'/®/47r0A2(0“^), where K 2 is the modified Bessel function of the 
second kind of order 2. The moments of the distribution are written as f Anp'^ f{p)dp = 1, 


J ^T^p‘^f{p)dp 

[ 4:Trp^f{p)dp 


30 + f0^ 


150^ 
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This correction arises another second-order spontaneous term on the order of A5(p4) in 
addition to Eq. (16) because we have 
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where the second term of the last expression is on the order of A5(p4) while the first term is 
already appeared in Eq. (12). 

Combining all of them, we obtain the higher-order Kompaneets equation that satisfies 
conservation of photon number, and the Bose-Einstein distribution is the equilibrium solution 
for this equation (c.f., the study for Eq. (12) is given by Caflisch and Levermore (1986) [33]). 
Note that the resultant equation is exactly the same as Eqs. (15) and (16) of Challinor and 
Lasenby (1998) [34], although they took a different approach to derive the equations. 


3. The Case for Large Occupation Number: Steady-State Solution and Its 
Stability 

We consider the case for large occupation number n ^ 1. When we start from the first order 
(Eq. (12)), for re S> 1, we obtain 

= 0 . ( 22 ) 

where N{x) = x^n{x). This is the nonlinear convection equation, which leads to an unphys¬ 
ical multi-valued solution. The situation is completely analogous to the shock formation in 
hydrodynamics so that the presence of a viscosity (the diffusion term in the right-hand side) 
can be the resolution (c.f. Ref. [35]). The viscous term arises when we extend the ideal fluid 
(the Euler equations) to the viscous fluid (the Navier-Stokes equations) (c.f. Ref. [36]). On 
the other hand, in this case of ICS, we will see that the physical resolution is not a diffusion 
term. 

Now, we start from the higher-order Kompaneets equation and find a new equation that 
resolves the difficulty of Eq. (22). Because we expand Eq. (1) assuming /c <C 1 and 0 <C 
1 (p <C 1), we neglect the second-order spontaneous terms (Eqs. (15) — (17) and (21)) 
compared with the first-order spontaneous terms and we obtain 
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Efficiencies of each term on the right-hand side of Eq. (23) are characterized by scattering 
optical depths r. In order of magnitude calculation, dividing each term by n/y, the right-hand 
side of Eq. (23) includes four components whose scattering optical depths are described as 
Tspon = xy (the first term: spontaneous Compton), = y (the second term: inverse Comp¬ 
ton), Tjjjj = xyn (the third term: induced Compton) and Tggc = xynQ <C other 

terms proportional to 0: second-order induced Compton), respectively. ICS becomes dom¬ 
inant process if ^ '^spoinRnvJ ^ ^ min(l,x) S> 1. If we consider only the term, 
the mathematically troublesome Eq. (22) is derived. On the other hand, when the 
terms are neglected, we obtain the usual Kompaneets equation Eq. (12). For the case of 
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X = h V / kg Tpi 

Fig. 1 Plots of Eq. (26) in logarithmic scale for different parameters: {kQ,A,B,(j)) = 
(tt, 5,5,0) (red line) and (dvr, 45, 55, 0) (blue line). The plasma temperature of the blue line 
is colder than that of the red one. 


Tgec/max(rspon) "^inv) ~ ^ mm(l,x) ^ 1, the first-order spontaneous terms are neglected 

compared with the second-order induced terms even though the terms are higher-order 
in A: or 0. This is because the spontaneous and induced terms have different dependence on 
n so that this does not mean an importance of the third and further higher-order terms in 
k or 0. 

The large occupation number satisfying the condition n0 x min(l,a:) = mm(nQ,nk) ^ 
1 is expected in some situations, for example, n ~ 10^^ (brightness temperature Tb(^) = 
hi'n{u)/kB ~ lO^^K at frequency u ~ lOGHz) is implied from an observation of the Crab 
pulsar in astrophysics [37] and n ~ 10^^ (petawatt laser with ~ 3.3 nm spectral width at 
wavelength ~ 10.5/rm) is obtained from some current laser facilities (e.g. Ref. [38]). In this 
case, we find that the natural extension of Eq. (22) is written as 


dN 


d 


14 


dy 


dfi 


17,, d 


— ^ I - -,N—N 


dx^ 


14 dx 


(24) 


Note that Eq. (24) exactly corresponds to the isotropic case of Eq. (4) with A5 = 0 and 
A/ = ^I{k) + contrast to Eq. (22), the plasma temperature 0 appears explicitly 

on the right-hand side that corresponds to A/(fcp 2 ) (Eq. (18)). Because the right-hand side is 
the higher-order correction term, the velocity of photon flux (in frequency-space) is basically 
determined by the second term of the left-hand side, i.e., order N{x) to the negative direction 
of X. The first term of the right-hand side is important to resolve the problem described in 
Section 1, while the second term slightly modifies the second term of the left-hand side. Eq. 
(24) is similar to the Korteweg-de Vries equation, but the coefficient of the dispersive term 
is also non-linear (proportional to N). This is different from the study by Peyraud (1968) 
[ 20 - 22 ]. 
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Before proceeding numerical calculations, let us find the steady-state solution. The exact 
analytic solution of Eq. (24) is obtained in steady state N(x), i.e., 
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(25) 


where we introduce the Doppler wavenumber fc© = (5/(70) — 31/14)^/^ ^ 1 for 0 <C 1. We 
obtain a non-trivial solution [N ^ 0) from the first line of Eq. (25) as 


N{x) = A cos{kQ In X -)- (/) + B, 


(26) 


where an amplitude A{> 0), a DC component B and a phase (j) are constants of integration. 
Eq. (26) shows no characteristic photon energy x, because the first-order spontaneous terms 
are neglected in Eq. (24), where only has different dependence on x from rspon,rjj^j 
and Tgg^. This steady-state solution N{x) includes the solution of Eq. (22), as A = 0, i.e., 
the first term of Eq. (26) is the important contribution from We require B > A 

for N(x) > 0. Inside the derivative of the second line of Eq. (25) is written as a constant 
(70/5)/c0(i?^ — A^) > 0, i.e., photon flux in frequency-space is negative direction of x. The 
photon flux is zero for A = B, i.e., contribution from compensates ^l(k) term in this 

case and plays as a kind of heating, while Al^k) term plays a role of cooling (see Eq. (13)). 
We require proper boundary conditions at finite frequencies to determine the constants of 
integration, because the number density f N{x)dx and the energy density f xN{x)dx of 
photons diverse for either x —> 0 or x —>■ 00. Since the photon flux is the negative direction, 
a photon injection at a large x and absorption at a small x are expected to achieve the 
steady-state (c.f., a similar study for Eq. (12) is given by Dubinov (2009) [39]). We plot 
Eq. (26) in Eig. 1, for example. 0ne of important features of Eq. (26) is that the Doppler 
wavenumber A:© ~ 0“^/^ associates with logarithmic scale in frequency x. This arises from 
the third-derivative term in Eq. (24) and is interpreted as the Doppler width Av rsj 01/2 

V 

by induced-inverse Compton process (A/(fcp 2 )) already discussed by Zel’dovich and Sunyaev 
(1972) [25] in the integral form. Note that this Doppler width is a different process from the 
first-order spontaneous inverse Compton scattering AS'(p 2 ). For the later convenience, we 
define the Doppler wavelength A 0 = 27r/A:©. For every integer m, the peaks of N{x) appear 
at X = with a constant <I> = 

Finally, we analyze linear stability for the steady-state solution N{x). We substitute 
N{x,y) = N{x) + Ni{x,y) {Ni/N <C 1) into Eq. (24) and linearize the equation. 


dNi 

dy 


14 / 93 


—QN 
5 




e + l) 


dNi 

dx 


(27) 


Considering a Fourier component of Ni{x,y) oc we obtain the dispersion relation 

oj = —QN {ixK^ — K^x^ — (k 0 -|- 1)k) . (28) 

5 

Because Im(a;) > 0, the steady-state solution Eq. (26) is linearly stable. 
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4. Numerical Simulation 

Here, we numerically study the initial evolution of photon spectra by solving Eq. (24) with 
n ^ 1. We pay particular attention to evaluate the third-derivative in Eq. (24) precisely in 
order to see developments of solitary structures quoted by Zel’dovich and colleagues [25, 26]. 

4-1. Set up 

Eq. (24) is solved by fourth-order in x derivative (five-point stencil) and fourth-order in y 
derivative (fourth-order Runge-Kutta). We use x grid points spaced linearly from 0.1 to 10 
and put fixed boundary conditions at both sides of x boundaries. Initial spectra are assumed 
to be a Gaussian spectrum, 

No{x) = Wnitexp (29) 

V / 

including three parameters, a normalization Njnit, a mean Xinit and a half width Uinit of a 
Gaussian. The ratio of the initial spectral width cjinit to the Doppler wavelength A0 = 27r/A:©, 
which is a parameter associated with Eq. (24), characterizes the spectral evolution. While suf¬ 
ficiently broad spectra (uinit A 0 ) are expected in astrophysical situations, spectral widths 
of laser beams in laboratories are not always broad. We choose initial conditions that satisfy 
^T-init© X min(l,Xinit) 3> 1 (xfjjj^ninit = Ninit). Since the optical depth to ICS is written as 
Tnd = yN/x, we measure time y with yo = a;init/Ninit and an adequate time-step is searched 
to resolve spectral evolution for each calculation. 

Eq. (24) is approximated form of Eq. (23) for n0 x min(l, x) ^ 1. We checked that both of 
the equations give almost the same initial evolution until y < 0.4yo for the same parameter 
sets, satisfying n0 x min(l, x) ~ 10 and 10^ at the peak of a Gaussian x = Xjnit- For 0.4yo ^ 
y < 2 / 0 ) numerically unstable structures start to develop at the portion of N{x,y) <C Ni^it 
in our current scheme. This does not allow us to pursue the spectral evolution for different 
(Tinit and A 0 beyond y r\j 0.42/0- In what follows, we show numerical results of Eq. (24) for 
y < OAyo. We fixed (iVinit, Xinit) = (10^, 1), since different sets of (Njnit, ®init) does not change 
results as long as we measure y with 2 / 0 - For remaining two parameters, we study the case 
that both (Tinit and A 0 are an order of 10“^ (A 0 = 0.1 corresponds to an electron temperature 
of ~ 10^ eV). 

Eq. (24) has two invariants of motion, the photon number density f N{x,y)dx and a 
quantity f \nN(x, y)dx [25]. Conservation of photon number is immediately found when we 
rewrite x derivatives in Eq. (24) to the second line expression of Eq. (25) and integrate 
both sides the equation over dx. We obtain conservation of the quantity f lnN(x, y)dx from 
dividing both sides of Eq. (24) by N{x) and integrating over dx. For all the results of our 
calculation below, we checked that these quantities were conserved for y < 0.42/0- 

4-2. Results 

Fig. 2 shows results of the numerical simulations for different Doppler widths A 0 = 0.20,0.10, 
and 0.05 where the initial spectral widths are (Tinit = 0.2 for the left panel and (Tinit = 0.1 
for the right panel, respectively. Black lines in Fig. 2 show the initial spectra Nq{x), and 
red and blue lines are photon spectra N{x,y) at y = 0.27/0 and at y = 0.47/0, respectively. 
We Hnd that solitary structures are formed intermittently and shift to lower frequency, i.e., 
direction of their motion is basically determined by the left-hand side of Eq. (24) and their 
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Fig. 2 Each panel shows results of numerical calculation of Eq. (24) for three different 
Doppler width A 0 = 27r/A:0 = 0.20 (uppermost lines), 0.10 (middle lines), and 0.05 (low¬ 
ermost lines). The initial spectral widths are different between left (iTinit = 0.2) and right 
panels (iTinit = 0.1), while (IVinit, a^init) = (lO'^jl) is common (see Eq. (29)), i.e., yo = 10“'^ 
for all cases. Black, red and blue lines correspond to A^(x,0), N{x,0.2yo) and N{x,OAyo), 
respectively. Insets are zoom of solitary structures appeared in the middle line (A© = 0.10) 
and the lowermost line (A© = 0.05) on the left panel and in the lowermost line (A© = 0.05) 
on the right panel when y = OAyQ. Blue lines in the insets are the result of calculation 
A^(x,0.4?/o) and green lines are Eq. (30) with parameters tabulated in Table 1. Numbering 
of solitary structures corresponds to the number i on Table 1. 


Table 1 Fitted values of rii and rrii in Eq. (30). Number i of solitary structures is found in 
the insets in Fig. 2. Corresponding dinit and A 0 are also tabulated. Arrii is separation between 
the neighboring peaks of the i-th solitary structure, Amia = uT, 2 a — for example. 


i 

finit 

II 

m 

rrii 

Arrii 

la 



3.42 X 10^ 

-0.869 

0.192 

2a 

0.20 

0.10 

2.38 X 10^ 

1.38 X 10^ 

-0.677 

0.173 

3a 

-0.504 

0.153 

4a 



4.78 X 10® 

-0.351 

- 

lb 



3.64 X 10 '^ 

-0.921 

0.0996 

2b 



3.14 X 10 '^ 

-0.821 

0.0960 

3b 



2.64 X 10 '^ 

-0.725 

0.0915 

4b 

0.20 

0.05 

2.16 X 10^ 

-0.634 

0.0873 

5b 



1.68 X 10^ 

-0.547 

0.0838 

6b 



1.15 X 10^ 

-0.463 

0.0795 

7b 



6.27 X 10® 

-0.383 

- 

Ic 



3.30 X 10'^ 

-0.481 

0.0981 

2c 

0.10 

0.05 

2.44 X 10^ 
1.68 X 10 '^ 

-0.383 

0.0917 

3c 

-0.292 

0.0844 

4c 



9.40 X 10® 

-0.207 

— 
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solitary form results from the first term of the right-hand side of Eq. (24). The heights of 
each solitary structure increase with time without changing their logarithmic width. The 
smaller value of A 0 , i.e., the colder plasma, the more and the narrower solitary structures 
are formed. This behavior is consistent with what is inferred from the steady-state solution 
(Eq. (26)). On the other hand, numerical simulations without the third-derivative in Eq. (24) 
gives an unstable shock-like discontinuous structure rather than solitary structures with the 
use of the same scheme. 

To characterize the solitary structures, we consider a zeroth-order log-normal distribution 
for each solitary structure and fit some of the results at y = 0.4yo with a function 

F{x) = Noix) + Y^ni exp ^ , (30) 

where we adopt a half width of u; = Ag/d considering Eq. (26). Eitted values of rii, rrii 
and the separation between the neighboring peaks Arm (for example, Amia = uT, 2 a — mia) 
are listed in Table 1. Numbering of solitary structures is found on the insets in Eig. 2; 
i = la — 4a for (uinit, Aq) = (0.20,0.10), i = lb — 7b for = (0.20,0.05) and i = Ic — 4c for 
= (0.10,0.05), respectively. We do not fit all the structures because some structures have 
a smaller amplitude than Nq{x). Eor example, we see more structures on the right of the 
structure numbered ‘4c’ of the inset on the right panel in Eig. 2. Logarithmic widths of 
the structures are well characterized by the Doppler width A 0 . However, the spacing of 
neighboring peaks Am, does not follow Eq. (26), i.e.. Ami 7 ^ A 0 (see Table 1), or rather 
Ami is larger for lower frequency peaks. This relates with larger m at lower frequency, 
because the velocity of a wave is proportional to N{x) in Eq. (24). Eitted values mia < mm 
indicate that the velocity of solitary structures is slower for larger value of A 0 . Although 
the velocity relates with energy loss of photons, it is not easy to discuss the dependence on 
(o'init) Ae) with Eig. 2, because number and height of structures are also different for different 
sets of (o-init,A 0 ) (see the discussion about photon energy transfer in Section 4.3). 

4.3. Discussion 

Radio emissions from pulsars (e.g. Ref. [37]) and also fast radio bursts (e.g. Refs. [41]) have 
extremely large brightness temperature. Although the optical depth to ICS rind would attain 
to unity for these objects (c.f. Refs. [15, 16, 18]), the past studies did not discuss what kind 
of signatures are expected to be imprinted in their observed spectra. Our present study 
predicts a spectral break at the frequency t'ind corresponding to rind(Dnd) ^ 0.2 and solitary 
structures below ^nd in tke photon spectra (red lines in Eig. 2). In this regard, the discrete 
emission bands in the dynamic spectra of the giant radio pulse occurred at the interpulse 
phase in the Crab pulsar reported by Hankins & Eilek [37] (‘zebra bands’) are intriguing 
phenomena. If we interpret their reported value Alnz/ = Av/v ~ 0.06 by ICS, the value is 
not far from Am* for Aq = 0.05 in Table 1 corresponding to the electron temperature of 
~ a few xlO eV. However, for applications to realistic astrophysical situations, we should 
take account for anisotropic photon distributions and for relativistic effects of both bulk and 
thermal motions of plasmas. These effects will be studied in subsequent papers. 

While the quantities f N{x,y)dx and J In N{x,y)dx are conserved with time, the photon 
energy density eph(y) = f xN{x,y)dx decreases with time, i.e., the energy of photons is 
transfered to plasmas by ICS [32, 40]. The rate for the energy transfer by ICS is estimated 
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y/vo 


Fig. 3 Evolution of the normalized photon energy density eph(?/)/eph(0) with the nor¬ 
malized time y/yo- The vertical axis is the natural logarithm of eph(y)/eph(0). Thick and 
thin lines correspond to (Tinit = 0.20 (left panel in Fig. 2) and = 0.10 (right panel in Fig. 2). 
We also show the results for cjinit = 0.05 (dotted lines) for two different Doppler widths of 
A 0 = 0.10 (red) and 0.05 (blue), respectively. 

from the integration of Eq. (22), dep^{y)/dy = — f N‘^{x,y)dx (see Eq. (13)). Considering 
the initial phase of evolution y ^ yo, i.e., N{x,y) ks Nq{x), we estimate the rate as 

c^eph(y) ,,,2 

~ for y < yo- (31) 

yo 

where we approximate eph(y) = / xN{x,y)dx « 2xinitO'init.A^init and recall yo = a^init/-A^init for 
the last expression. Eq. (31) can be applicable to describe an initial evolutionary phase and 
gives exponential loss of the photon energy with time eph(y) ~ eph(O) exp(—y/yo). 

Fig. 3 plots evolution of the normalized photon energy density with time for the results 
of Fig. 2. We also show the case for an initial width of cJinit = 0.05 as dotted lines in Fig. 3. 
Thin lines (dinit = 0.20) show the exponential energy loss, although the slope is not exactly 
ln[eph(y)/eph(0)] = —y/yo- Thick (dinit = 0.10) and dotted (dinit = 0.05) lines deviate from 
the exponential-law (thin lines) at y ~ 0.15yo and y ~ 0.05yo, respectively. Dotted red line 
(TnitjAe) = (0.05,0.10) and thick black line (dinit, Ae) = (0.10,0.20) show a bit different 
slope even at y < 0.05yo from other lines and these exceptional behaviors are characterized 
by a relatively large Doppler wavelength Ae = 2dinit, which will be discussed later in this 
section. The energy transfer from photons to plasmas by ICS is effective for the case of broad 
initial spectra than narrow ones (the blue lines are always below the red lines for the same 
o'init)- This means that the deviation of N{x,y) from the initial spectrum No(x) is faster 
for narrower spectra and this behavior is found from the inset on the right panel of Fig. 
2 (the blue line is well below the green line compared with those on the left panel of Fig. 
2). Although we see a temperature dependence of the energy transfer for the same initial 
width, Eq. (31) has no information of electron temperature explicitly. The third-derivative 
in Eq. (24) may play as a heating term because larger Ae (higher temperature) gives slower 
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energy transfer in Fig. 3. This is consistent with the discussion at the steady-state solution 
Eq. (26), i.e., the photon flux is written as ~ for 0 <C 1 and a contribution from 

A/(fcp 2 ) (A) resists that from (B). 

We should note the case of Ae = 2cJinit, i-e., full width of the initial spectrum is almost the 
same as the Doppler width. For Ae > 2(Tinit, evolution of photon spectra is unstable at least 
in our numerical scheme. In those calculations, although all of the results shown in this paper 
have no time variation of N(x > Xinit, y) (see Fig. 2), there appears unstable features at x > 
Xinit for A© > 2crinit- This unstable feature is not improved by higher numerical resolutions 
of both X and y. Currently, it is uncertain whether the exceptional behaviors of dotted 
red line ((Tinit,A 0 ) = (0.05,0.10) and thick black line (uinitjA©) = (0.10,0.20) in Fig. 3 is 
numerical or physical. However, we should also take care of our formulation (Eq. (1)) for a 
narrow spectrum Uinit ‘C A©. Galeev and Syunyaev (1973) [ 6 ] argued that, when Uinit ^ A©, 
collective behaviors of plasmas are important rather than the interaction with free electrons, 
i.e., Compton scattering. This is the reason why we do not study the case for Aq > 2iTinit more 
in this paper, although study of this regime is important for the application to laboratory 
experiments. 

We should also note that whole of present paper is based on the assumption 0 = const. 
As seen in Fig. 3, photons clearly lose energy by ICS and transfer their energy to plasmas. 
The time-scale of energy transfer to electrons is estimated from Eq. (31) as 

f 1 depi{y) \~^ ^ / 1 dephiy) \~^ 

Vepi(y) dy J V epi(y) dy ) 

~ for y < yo, (32) 

eph(y) 

where Cpi ~ n^Q is the energy density of a plasma divided by mgC^. We need to take into 
account evolution of a plasma temperature for evolution beyond this time-scale (c.f. Refs. 
[23, 24]). Note that Eq. (32) is allowed to be used even for epi(y) > eph(y) (c.f. Ref. [7]). Cur 
calculation needs only Xinit, A^init, crinit and A© that determine yoi eph(O), and 0, i.e., we do 
not specify Ue and t separately. Typical time-scales of photon cooling tph (Eq. (31) and also 
Eq. (24)) is found from yph ~ yo and that of electron heating tpi (Eq. (32)) is found from 
ypi ~ yo(epi/eph)- These are expressed as 


tph ~ 


Xinit 1 , 1 TT Ag 

WnitC ngCTrpC’ 2crinitfVi^it0'^ O'rpC ’ 


(33) 


where the energy densities of photons and electrons are 6 *^^ = lic 0 ^/( 7 r^Ag) f xN{x)dx and 
€*[ ~ UeOmeC^ including dimension, respectively. For given Xinit, A^init, and A©, photons 
lose energy before heating electrons (tph fpi) for sufficiently large Ug. For example, we 
obtain tph < tpi 10 ps for Ug > lO^^cm ^ (close to conditions of laser experiments) with 
the parameters which we used in this paper. 

We finally discuss differences from past calculations of spectral evolution by ICS [28, 29]. 
Fig. 1 of Coppi et al. (1993) [29] is spectral evolution for isotropic case. Basically, they solved 
Eq. (22) but including numerical viscosity, i.e., their calculations did not have information 
of a plasma temperature (A©). Although their result (the right panel of Fig. 1 of Coppi et al. 
(1993) [29]) showed two solitary structures that have the same logarithmic width but their 
behaviors seem different from ours, such as the height of the low frequency structure is smaller 
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than that of the high frequency one. We consider that these structures are numerical, i.e., 
combination of the numerical viscosity and logarithmically-spaced frequency grids in their 
calculations. Montes (1979) [28] solved the integro-differential equation given by Zel’dovich 
et al. (1972) [26], but with some simplifications. Fig. 2 of Montes (1979) [28] (the calculation 
for an initially broad Gaussian spectrum) shows solitary structures of the linearly same 
width. We consider that Eq. (27) of Montes (1979) [28] is over simplified, especially, their 
assumption Au ~ = const, in their Eq. (20) is crucial, while we consider the case 

Av ^ [25]. 

5. Conclusions 

In this paper, we study evolution of photon spectra when ICS dominates. To get rid of the 
well-known difficulty of the nonlinear convection equation (Eq. (22)), we consider the higher- 
order Kompaneets equation. We obtain the new equation (Eq. (24)) that describes evolution 
of photon spectrum by ICS and that overcomes the difficulty. The second-order induced 
term (Eq. (18)) obtained from the higher-order Kompaneets equation improves the 

formulation. In addition, Eq. (24) has the steady-state analytic solution (Eq. (26) and Fig. 
1), which is linearly stable. The steady-state analytic solution predicts the formation of 
solitary structures of logarithmically same width in frequency-space and the A/(fcp 2 ) term 
plays as a heating term against the ICS term A/^^), which is the pure cooling term. 

We also study evolution of photon spectra by ICS numerically. ICS intermittently forms 
solitary structures moving toward lower frequency (Fig. 2) and these behaviors are con¬ 
sistent with predictions by some of past studies [25, 26]. The solitary structures have the 
logarithmically same width well characterized by the Doppler width A© and this behavior is 
also inferred from the steady-state solution. On the other hand, the spacing between peak 
frequencies of the structures does not follow Eq. (24) {Arrii ^ A©). The number and height 
of solitary structures depend on A©. 

The results of our numerical simulation satisfy the two invariants of motion, conservations 
of photon number f N{x)dx and of the quantity f In N{x)dx. On the other hand, the energy 
density of photons f xN(x}dx is transferred to electrons (Fig. 3). Energy density decays 
exponentially in an initial phase and this behavior is almost consistent with the analytic 
estimate (Eq. (31)). The energy transfer from photons to plasmas by ICS is effective for 
the case of broad initial spectra such as expected in astrophysical situations. The numerical 
results shown in the present paper are only initial phases of evolution y < OApo- We need 
more optimized numerical scheme to study evolution for y > 0.4yo- 
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